params <- matrix(nrow=1080,ncol=6)
dpl <- c(0,0.25,0.5,0.75,1)
scl <- c(0.001,0.01,0.1,0.5,1,2)
sml <- c(10^-5,10^-6,10^-7,10^-8)
nml <- c(10^-4,10^-5,10^-7)
sw <- 0

for (divprob in 1:5){
  for (selcoeff in 1:6){
    for (selmu in 1:4){
      for (neutmu in 1:3){
        for (replicates in 1:3){
          sw <- sw + 1
          params[sw,1] <- sw
          params[sw,2] <- dpl[divprob]
          params[sw,3] <- scl[selcoeff]
          params[sw,4] <- sml[selmu]
          params[sw,5] <- nml[neutmu]
          params[sw,6] <- replicates
        }
      }
    }	
  }
}

# Funky plot
# Set ylim on divergence
y.lim.low <- -0.20
y.lim.upp <- 0.05
col.scheme <- c("black", "blue", "darkgreen", "orange", "red")

pdf(paste("biopsy.sampling.pdf",sep=""),width=6,height=6,pointsize=12)
for (selcoeff in 1:6){
  for (selmu in 1:4){
    for (neutmu in 1:3){
      ind <- params[params[,3]==scl[selcoeff]&params[,4]==sml[selmu]&params[,5]==nml[neutmu]&params[,6]==1,1]
      error.div.est <- matrix(nrow=length(ind)*10,ncol=2)
      for (j in 1:length(ind)){
        i <- ind[j]
        divergence <- read.csv(paste("../sim",i,"/sim",i,".divergence.log",sep=""),header=FALSE)
        biopsy <- read.csv(paste("../sim",i,"/sim",i,".divergence.biopsy.sampling.log",sep=""),header=FALSE)

        real.div <- divergence[divergence[,1]==7300,2]
        for(biop in 1:10){

          f <- c(10,20)
          for(cellind in 1:2){
            cell <- f[cellind]
            div.est <- biopsy[biopsy[,2]==biop&biopsy[,3]==cell,]
            error.div.est[((j-1)*10)+biop,cellind] <- median(div.est[,4])-real.div
            #low.div.est[biop,cellind] <- as.numeric(quantile(div.est[,4])[1])
            #upp.div.est[biop,cellind] <- as.numeric(quantile(div.est[,4])[4])
          }
          
          
        }
      }
      plot(0,0,t="n",ylim=c(y.lim.low,y.lim.upp),xaxt="n",xlim=c(1,10),xlab="No. Biopsies",ylab="Divergence Estimate Error")
      title(main=paste("sc =",scl[selcoeff],"smu =",sml[selmu],"nmu =",nml[neutmu]))
      for (j in 1:length(ind)){
        lines(c(1:10),error.div.est[((j-1)*10)+1:10,1],lty=1,col=col.scheme[which(dpl==params[ind[j],2])])
        lines(c(1:10),error.div.est[((j-1)*10)+1:10,2],lty=2,col=col.scheme[which(dpl==params[ind[j],2])])
      }
      abline(h=0, col="black")
      axis(side=1,at=c(1:10),labels=c(1:10))
    }
  }
}
dev.off()
		
